#!/usr/bin/perl

# ------------------------------------------------------------------------------
# --- Use CD-HIT to obtain non-redundant entries 
# ---  Output is pdb_dssp_seq_95.fa 
# --- PDB codes of nonredundant entries from CD-HIT *.clstr output are written 
# ---  to nonredundant_pdb.txt 
# --- get_fasta_from_dssp.pl must be run first!
# ------------------------------------------------------------------------------

use strict;
use warnings;

# Be sure to run "module load perl" from the PBS script first.
use Bio::DB::Fasta;
use Bio::SeqIO;

print STDERR "Calling CD-HIT to obtain non-redundant entries...\n";

my $db = Bio::DB::Fasta->new('ss.txt');

my $cd_hit_path = "cd-hit-v4.5.4-2011-03-07/";
my $cd_hit_cmd = $cd_hit_path . "cd-hit -i pdb_dssp_seq.fa ";
$cd_hit_cmd .= "-o pdb_dssp_seq_95.fa -c 0.95 -n 5  > cd-hit.out 2>&1;";

system($cd_hit_cmd);

# This creates two files, a FASTA file of the non-redundant sequences 
#  (pdb_dssp_seq_95.fa) and info on the clusters (pdb_dssp_seq_95.fa.clstr).

my $seq90_in = Bio::DB::Fasta->new('pdb_dssp_seq_95.fa');
my @seq90_ids = $seq90_in->ids;

# Print out number of non-unique sequences
my $dssp_seq_count_nonuniq = `grep '>' pdb_dssp_seq_95.fa | wc -l`;
$dssp_seq_count_nonuniq =~ s/\s//g;

print STDERR "\tNumber of non-unique sequences found in DSSP file:\t$dssp_seq_count_nonuniq\n";

# ------------------------------------------------------------------------------
# --- Grab PDB codes of nonredundant entries from CD-HIT *.clstr output --------
# ------------------------------------------------------------------------------

my $grep_cmd = "grep '*' pdb_dssp_seq_95.fa.clstr | ";
$grep_cmd   .= "cut -d' ' -f 2 | cut -d':' -f1 | cut -d'>' -f2 | ";
$grep_cmd   .= "sort | uniq > nonredundant_pdb.txt";

system($grep_cmd);

# Print out number of unique sequences
my $dssp_seq_count_uniq = `wc -l nonredundant_pdb.txt`;
$dssp_seq_count_uniq =~ s/\D//g;

print STDERR "\tNumber of unique sequences found in DSSP file:\t$dssp_seq_count_uniq\n";
